# ------------------------------------------------------------------------------------------------
### Generate Figure 18: Effects for municipalities <300,000 inhabitants
# ------------------------------------------------------------------------------------------------

# Load data
load("../data/results/sdid_results.Rda")

# Define the outcomes to plot with their display names
outcome_mapping <- c(
  "Opposition" = "Opposition Turnout",
  "Law and Justice (PiS)" = "Government Turnout", 
  "Overall Turnout" = "Overall Turnout"
)

# Filter and prepare data 
plot_data <- sdid_results %>%
  filter(
    # Filter for main analysis (Any Level, Any Type, notyettreated, No cities >300k, Full Sample)
    sum_level == "Any Level" & 
    sum_type == "Any Type" & 
    level == "Any Level" & 
    type == "Any Type" & 
    control_sample == "notyettreated" & 
    no300k_sample == "No cities >300k" & 
    noproposed_sample == "Full Sample",
    # Filter for main outcomes
    (party == "Opposition" & outcome == "Party Turnout") |
    (party == "Law and Justice (PiS)" & outcome == "Party Turnout") |
    (party == "Overall Turnout" & outcome == "Overall Turnout")
  ) %>%
  mutate(
    outcome = outcome_mapping[party],
    outcome = factor(outcome, levels = c("Government Turnout", "Opposition Turnout", "Overall Turnout"))
  )

# Generate plot
ggplot(plot_data, aes(x = estimate, y = outcome, color = sample)) +
  # Confidence intervals (99%, 95%, 90%) with decreasing linewidth
  geom_errorbarh(aes(xmin = lower_cb_99, xmax = upper_cb_99), height = 0, linewidth = 1, alpha = 0.5, position = position_dodge(0.3)) +
  geom_errorbarh(aes(xmin = lower_cb_95, xmax = upper_cb_95), height = 0, linewidth = 3, alpha = 0.3, position = position_dodge(0.3)) +
  geom_errorbarh(aes(xmin = lower_cb_90, xmax = upper_cb_90), height = 0, linewidth = 5, position = position_dodge(0.3)) +
  # Point estimates
  geom_point(shape = 21, fill = "white", size = 3, stroke = 1.5, position = position_dodge(0.3)) +
  # Styling
  scale_color_manual(values = c("Full Sample" = "black", "< 50km Sample" = "grey"), guide = guide_legend(title = NULL)) +
  xlim(-2, 2) +
  labs(x = "Percentage points") +
  geom_vline(xintercept = 0, linetype = "dotted") +
  theme_bw() +
  theme(
    legend.position = "bottom", 
    axis.title.y = element_blank(), 
    panel.grid.major.y = element_blank(), 
    text = element_text(size = 23)
  )

# Save the plot
ggsave("../output/plots/figure18.pdf", width = 14, height = 7, dpi = 700)

